# Replication Script for BJPS-Article "Setting the party agenda: Interest groups, voters and issue attention"
# by Heike Kl�ver, 28.12.2017
# ============================================================================================================
  
# load libraries
library(foreign)
library(nlme)
library(lattice)
library(Matrix)
library(lme4)
library(car)
library(arm)
library(memisc)
library(mice)
library(Zelig)
library(Amelia)
library(plotrix)   
library(graphics)     
options(scipen=3)       
library(readstata13)



# =========================================
# Figure 1  (based on model 2)
# =========================================

# clear working space
rm(list=ls(all=TRUE)) 
# set working directory
setwd("C:/Users/")

# read in coefficients
fixeff_raw <- read.csv("fixeff_m1.csv")
# select first row
fixeff_untransposed <- fixeff_raw[1,]
# transpose vector
fixeff2 <- t(fixeff_untransposed)
# remove row names
fixeff <- fixeff2[1:12]
# read in variance covariance matrix
varcov_raw <- read.csv("varcov_m1.csv")
# retain the first 12 x 12 matrix
(varcov <- varcov_raw[1:12,1:12])



# Step 1: get b and V
# ====================
# create vector b with coefficient
(b <- as.vector(fixeff))		  
# create matrix v with variance-covariance matrix
(V <- as.matrix(varcov))		
# glimpse at the standard errors for each coefficient
(sqrt(diag(V)))                         



# Step 2: Set up MVN distribution
# ===============================
# choose number of simulations
nsim <- 1000
# draw from multivariate normal distribution
# with nsim number of simulations, 
# the coefficients b and variance-covariance matrix V
mvn <- mvrnorm(nsim, mu=b, Sigma=V)         


# Step 3: Get coefficients
# ========================
# read in data for figure 1 based on model 2
dat <- read.dta13("data_m1.dta")
attach(dat)
str(dat)


# Predicted values
# ================

# initialise empty list to be filled
ig_pred <- list(rep(NA,16))
# create vector containing sequence of numbers from 1 to 15
v_ig <- seq(0, 15, by=1)
# iterate through values 1 to 16
# to obtain linear predictions at
# different values of variable of interest
# and save results in empty list created above
for (i in seq_along(v_ig)){
  ig_pred[[i]] <-  (
    + mvn[,1] * i
    + mvn[,2] * mean(dat$lag_prev_mip_voter)
    + mvn[,3] * i * mean(dat$lag_prev_mip_voter)
    + mvn[,4] * mean(dat$lag_prev_finmin_rel)
    + mvn[,5] * mean(dat$lag_pervote)
    + mvn[,6] * mean(dat$extreme)
    + mvn[,7] * median(dat$gov)
    + mvn[,8] * mean(dat$lag_system_salience)
    + mvn[,9] * mean(dat$unemploy)
    + mvn[,10] * median(dat$unific)
    + mvn[,11] * mean(dat$lag_salience)
    + mvn[,12] * 1
  )
}



# Transform lists into matrixes
# =============================

ig_pred_matrix <- do.call(cbind, ig_pred)     



# Taking the means and quantiles of simulated values
# ==================================================

# initialise empty lists to be filled with quantities of interest below
ig_pred_mean <- list(rep(NA,16))
ig_pred_lb <- list(rep(NA,16))
ig_pred_ub <- list(rep(NA,16))

# obtain mean value and 95% confidence interval values 
for (i in seq_along(v_ig)){
  ig_pred_mean[[i]] <- mean(ig_pred_matrix[,i])
  ig_pred_lb[[i]] <- quantile(ig_pred_matrix[,i], c(.025))
  ig_pred_ub[[i]] <- quantile(ig_pred_matrix[,i], c(.975))
}




# Figure 1 (Interest groups)
# ==========================

win.metafile(file="C:/Users/fig1.emf", height=8, width=11)
par(mfrow=c(1,1),mar=c(5,5,1.5,3)) 
plot(1, type="n",  xlim=c(0,15), ylim=c(2,10),  xlab= expression("Interest group mobilization"~t[-1]),
     ylab=expression("Party issue attention"~t[0]~"(incl. 95% CI)"), cex=1.7, cex.axis=1.5, cex.lab=1.5)
lines(0:15,ig_pred_mean, col="black", lwd=2)
lines(0:15,ig_pred_ub, col="black", lty=2, lwd=2)
lines(0:15,ig_pred_lb, col="black", lty=2, lwd=2)
rug(dat$lag_prev_ig_rel)
dev.off()



# =========================================
# Figure 2 (Marginal effect plot)
# =========================================


rm(list=ls(all=TRUE))
setwd("C:/Users/")

# read in data for figure 2
dat <- read.dta13("figure2.dta")
attach(dat)
str(dat)

# read in data for figure 2 based on model 2
dat2 <- read.dta13("data_m1.dta") 
str(dat2)


win.metafile(file="C:/Users/fig2.emf", height=8, width=11)
par(mfrow=c(1,1),mar=c(5,5,1.5,3))
plot(1, type="n",  xlim=c(0,20), ylim=c(0,1),  xlab= expression("Party supporter issue attention"~t[-1]),
     ylab=expression("Marginal effect of interest group mobilization"~t[-1]~"(incl. 95% CI)"), cex=1.7, cex.axis=1.5, cex.lab=1.5)
lines(0:20,conbx, col="black", lwd=2)
lines(0:20,upper, col="black", lty=2, lwd=2)
lines(0:20,lower, col="black", lty=2, lwd=2)
rug(dat2$lag_prev_mip_voter)
dev.off()






# =========================================
# Table 2  (based on model 2)
# =========================================

rm(list=ls(all=TRUE))
setwd("C:/Users/")

# read in coefficients
fixeff_raw <- read.csv("fixeff_m1.csv")
# select first row
fixeff_untransposed <- fixeff_raw[1,]
# transpose vector
fixeff2 <- t(fixeff_untransposed)
# remove row names
fixeff <- fixeff2[1:12]
# read in variance covariance matrix
varcov_raw <- read.csv("varcov_m1.csv")
# retain the first 12 x 12 matrix
(varcov <- varcov_raw[1:12,1:12])


# Step 1: get b and V
# ====================
# create vector b with coefficient
(b <- as.vector(fixeff))		  
# create matrix v with variance-covariance matrix
(V <- as.matrix(varcov))		
# glimpse at the standard errors for each coefficient
(sqrt(diag(V)))     


# Step 2: Set up MVN distribution
# ===============================
# choose number of simulations
nsim <- 1000
# draw from multivariate normal distribution
# with nsim number of simulations, 
# the coefficients b and variance-covariance matrix V
mvn <- mvrnorm(nsim, mu=b, Sigma=V)       

# Step 3: Get coefficients
# ========================

dat <- read.dta13("data_m1.dta")
attach(dat)
str(dat)


# Predicted values --- change lag_prev_mip_voter from minimum, mean to maximum value to get the different predicted values and rerun again
# ================
# initialise empty list to be filled
ig_pred <- list(rep(NA,16))
# create vector containing sequence of numbers from 1 to 15
v_ig <- seq(0, 15, by=1)
# iterate through values 1 to 16
# to obtain linear predictions at
# different values of variable of interest
# and save results in empty list created above
for (i in seq_along(v_ig)){
  ig_pred[[i]] <-  (
    + mvn[,1] * i
    + mvn[,2] * mean(dat$lag_prev_mip_voter)
    + mvn[,3] * i * mean(dat$lag_prev_mip_voter)
    + mvn[,4] * mean(dat$lag_prev_finmin_rel)
    + mvn[,5] * mean(dat$lag_pervote)
    + mvn[,6] * mean(dat$extreme)
    + mvn[,7] * median(dat$gov)
    + mvn[,8] * mean(dat$lag_system_salience)
    + mvn[,9] * mean(dat$unemploy)
    + mvn[,10] * median(dat$unific)
    + mvn[,11] * mean(dat$lag_salience)
    + mvn[,12] * 1
  )
}



# Transform lists into matrixes
# =============================

ig_pred_matrix <- do.call(cbind, ig_pred)     



# Taking the means and quantiles of simulated values
# ==================================================
# initialise empty lists to be filled with quantities of interest below
ig_pred_mean <- list(rep(NA,16))
ig_pred_lb <- list(rep(NA,16))
ig_pred_ub <- list(rep(NA,16))


# obtain mean value and 95% confidence interval values 
for (i in seq_along(v_ig)){
  ig_pred_mean[[i]] <- mean(ig_pred_matrix[,i])
  ig_pred_lb[[i]] <- quantile(ig_pred_matrix[,i], c(.025))
  ig_pred_ub[[i]] <- quantile(ig_pred_matrix[,i], c(.975))
}





# ===================
# Predicted values
# ===================

# 1 (0) - due to the indicator variable, the first vector corresponds to an ig mobilization value of 0
mean(ig_pred_matrix[,1])
quantile(ig_pred_matrix[,1], c(.025, .975))

# 6 (5)
mean(ig_pred_matrix[,6])
quantile(ig_pred_matrix[,6], c(.025, .975))

# 11 (10)
mean(ig_pred_matrix[,11])
quantile(ig_pred_matrix[,11], c(.025, .975))

# 16 (15)
mean(ig_pred_matrix[,16])
quantile(ig_pred_matrix[,16], c(.025, .975))

